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Abstract We review the recently proposed extension of the Gutzwiller approxima- 
tion, M. Schiro and M. Fabrizio, Phys. Rev. Lett. 105, 076401 (2010), designed 
to describe the out-of-equilibrium time-evolution of a Gutzwiller-type variational 
wave function for correlated electrons. The method, which is strictly variational 
in the limit of infinite lattice-coordination, is quite general and flexible, and it is 
applicable to generic non-equilibrium conditions, even far beyond the linear re- 
sponse regime. As an application, we discuss the quench dynamics of a single-band 
Hubbard model at half-filling, where the method predicts a dynamical phase tran- 
sition above a critical quench that resembles the sharp crossover observed by time- 
dependent dynamical mean field theory. We next show that one can actually define 
in some cases a multi-configurational wave function combination of a whole set of 
mutually orthogonal Gutzwiller wave functions. The Hamiltonian projected in that 
subspace can be exactly evaluated and is equivalent to a model of auxiliary spins 
coupled to non-interacting electrons, closely related to the slave-spin theories for 
correlated electron models. The Gutzwiller approximation turns out to be nothing 
but the mean-field approximation applied to that spin-fermion model, which dis- 
plays, for any number of bands and integer fillings, a spontaneous Z? symmetry 
breaking that can be identified as the Mott insulator-to-metal transition. 



1 Introduction 

Time-resolved spectroscopies are advancing incredibly fast towards accessing ultra- 
short time (< femtoseconds) dynamics. HI |2] [3] |4]|5| On such timescales, it becomes 
possible to monitor how the electronic degrees of freedom react to a sudden external 
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Stimulus before electrons have time to equilibrate with the lattice, which commonly 
starts after few picoseconds. In this initial transient regime, one can therefore ne- 
glect the coupling to the lattice and study just the way how collisions among the 
electrons brought about by interaction redistribute the excess energy injected into 
the system. This situation in which the electrons provide their own dissipative bath 
has recently attracted interest especially in connection with cold atoms trapped in 
optical lattices, |j6l which realize systems where the particles are, to a large extent, 
ideally isolated from the environment. There are by now several claims that, when 
correlation is strong enough and the injected energy exceeds a threshold, the elec- 
trons alone are unable to efficiently exchange energy by collisions, hence remain 
trapped for long time into non-thermal configurations. The most convincing evi- 
dences come from dynamical mean field theory (DMFT) simulations of quantum 
quenches in the half-filled single-band Hubbard model. ||7] |8l Such a technique is 
however computationally heavy and does not allow accessing very long times. Al- 
ternatively, qualitatively similar results have been reproduced by a much simpler 
tool, the time-dependent Gutzwiller approximation (t-GA),||9l[l0l which allows to 
follow much longer the time evolution, although it lacks enough dissipative chan- 
nels to describe the system flowing towards a steady state. ifTol Nevertheless, the 
time averages of the observables as obtained through t-GA agree satisfactorily with 
the DMFT steady state values, which justifies using t-GA as a valid alternative to 
more sophisticated approaches, like DMFT, for its simplicity and flexibility. 

Here, we shall present in detail how t-GA can be implemented efficiently in a 
generic multi-band lattice model of electrons mutually coupled by a short-range in- 
teraction. We will show that the method is able to access the full out-of-equilibrium 
dynamics also far beyond the linear response regime discussed in Ref. ifTTI . In par- 
ticular, a nice feature of t-GA is its ability of treating on equal footing the dynamics 
both of the low-energy coherent quasiparticles as well as of the high-energy inco- 
herent excitations, which are commonly refereed to as the Hubbard side-bands close 
to the Mott transition. Within t-GA these two distinct excitations, quasiparticles and 
Hubbard bands, possess their own dynamics, and influence each other only in a 
mean-field like fashion. This is clearly an approximation of the actual time evolu- 
tion, and the reason why the method lacks enough dissipation, although the ensuing 
dynamics is much richer than the conventional time-dependent Hartree-Fock. 

Finally, we discuss some instructive connections between t-GA and the recently 
developed slave-spin representations of the Hubbard model. lfT2l [Tsl 1141 [Tsl Es- 
sentially, we will show that in the limit of infinite lattice-coordination, where the 
Gutzwiller approximation becomes an exact variational approach, and under partic- 
ular circumstances, e.g. integer filling in a multi-band model, one can actually define 
a multi-configurational basis of Gutzwiller wave-functions and explicitly evaluate 
the Hamiltonian matrix elements. It turns out that the Hamiltonian projected onto 
that basis coincides with its slave-spin representation with the major advantage that 
the constraint required in the slave-spin theory to project the enlarged Hilbert space 
onto the physical one can be here enforced exactly. 



The Out-of-Equilibrium Time-Dependent Gutzwiller Approximation 



3 



2 The model and the Gutzwiller wavefunction and 
approximation 

We shall consider the following tight-binding model on a lattice with coordination 
number z: 

t {!fj4.cj,+ii.c)+j:^, (1) 

i.j a.b=l i 

where c]^ creates an electron at site ; in orbital a = l,...,N, the index a includ- 
ing also the spin, and is a local term that accounts also for the interaction. The 
hopping parameter tfj' is assumed to scale like 1 /z' /^ where r is the lattice distance 
between sites ; and j, so that the average hopping energy per site remains finite also 
in the limit z — ^ °°. ||T6l The Gutzwiller wavefunction ifTTl [TSl is defined through 

\'F)^^\%)^]J^,\%), (2) 

i 

where | ff)) is a Slater determinan{3 and a local operator that we will denote, 
although improperly, as the Gutzwiller projector, whose role is to the change the 
weights of the local electronic configurations with respect to the Slater determinant. 
Both I %) and have to be determined variationally to minimize the total energy 

^ (f I 1 1') 



The Guzwiller approximation begins by imposing, for reasons that will become 
clear soon, the following two constraints on ^,: |[T9l 

(fol^/^,. |fo) = 1, (4) 
{% I 4,c., I %) = (fo I cLc,, I fo), Wa,b. (5) 

These constraints mean that, if we select from the operator any two fermionic 

operators and average over the Slater determinant what remains, then such an aver- 
age vanishes identically. This property is very convenient if the lattice coordination 
Z tends to infinity. In fact, we note that, for / 7^ j, 

{% I ^; ^, I fo) = {% I ^1^. I %){% I I fo) 

+ (f0 I ^l^i I fo)co„nected 

= l + {%\ I %>connected, (6) 

where the last term on the right hand side includes all Wick's contractions connect- 
ing the two sites, and the constant 1 comes from (|4|i. Because of the constraint (|5]l. 



^ In reality, for the method to work it is enough that Wick's theorem appUes, hence | ffj) could even 
be a BCS wavefunction. Here, for sake of simplicity, we shall only consider Slater determinants. 
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the terms that connect the two sites by only two fermionic Unes vanish, leaving only 
terms with 2« > 2 connecting lines. In the limit of infinite lattice-coordination, these 
latter terms vanish like z^"^'' , where Rij is the minimum length of the path connect- 
ing i to j. For a given /, if we consider all sites j at fixed Rij ~ R and sum over 
them Eq. (|6]l, each connected term above will contribute ^ z^"^, n > I, but there 
are only ^ z'^ such terms so that, in the limit z ^ °°, their sum will vanish. This 
property simplifies considerably all calculations in the infinite lattice-coordination 
limit, which we shall assume hereafter In particular, it implies that lfT9l l20l 

{^\W)=ll{%\^l^, |fo) = l, 

namely the wavefunction (|2]i is normalized, and moreover that, given any local op- 
erator ^i, 

(f |^Hf> = (fol^'^/^/ |fo), (7) 
which can be easily evaluated by Wick's theorem. In addition, it also follows that 

'■j <j 

where one has to keep only Wick's contractions that connect sites / and j by just 
a single fermionic line, since the terms with three or more lines vanish in the limit 
Z — >■ oo. A simple way to proceed is by defining the matrix elements Ri^b through 

{% I ^14,^. c,, I fo) ^ E Rl, {% I cj,c. , I fo), (9) 

c 

that automatically include all Wick's contractions after extracting from the operator 
^J4^i a single fermionic line. Through (|9]l we can formally write Eq. (O as 

('f' I ^i'^jh RiaRjbci ('f'o I clc., I . (10) 

cd 

In conclusion, provided (HI and Q are satisfied, and upon defining through 
Eq. ^ the renormalized hopping amplitude 

4^E4cf/;'^;.., (11) 

cd 

and the non-interacting Hamiltonian 

^*=EEO:t-^Ls-.+^-^-)' ^^^^ 

i,i ab 

then the average energy in the limit of infinite lattice coordination is 

£ = (fo I I fo> + E ('^'o I ^- \ %), (13) 
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which can be evaluated by Wick's theorem. Minimization of ([13) with respect to 
all variational parameters provides an estimate of the ground state energy. The ex- 
pression (flJl l, with the definition (|9]l, is strictly valid only in the limit of infinite 
lattice-coordination. However, it is common to keep using the same expressions 
also for finite-coordination lattices, hence the name Gutzwiller approximation. 

Like any other variational approach, also the one we just outlined can only pro- 
vide information on static properties, assumed to represent well those of the actual 
ground state. Here we shall propose an extension that allows to access also dynam- 
ical propertie s . Il9l [TOl 



3 Time-dependent Gutzwiller approximation 

From now on we shall assume that both the Slater determinant as well as the 
Gutzwiller projectors are time-dependent, hence 

I ^(t)) - 1 =n^'-(o I m))- (14) 

If the Eqs. (HI and Q are satisfied at any t, then, at any instant of time and in the 
limit of infinite coordination number, the average value of the Hamiltonian ) will 
have the same expression as in Eq. ( fTSl l. i.e. 

E{t) = {%it) I ^,(0 I %it)) I ^,-(f)t^-^,(0 I %{())■ (15) 

In particular, Mi{t) becomes time dependent since Riub{t) depends on time. We 
shall adopt the variational principle that | f (f)) is as close as possible to the 
solution of the Schroedinger equation. Specifically, [[9] we define the functional 
^{t) = JqcIt ^{t), that plays the role of a classical action, with Lagrangian 

J^it) = /(f^(f) I 'f'(f)) -£(f) = i{%{t) I ^(f)t^(f) I %{t)) 

+i{%{t) I J^(f)^^(0 I %{t))-E{t), (16) 

and determine | and ^,(f) by the saddle point of the action under the two 

constraints Eqs. (|4]i and Q- 

Since | is a Slater determinant at any instant of time, then 

^|%(0> = r(f)|fo(0), 

with 

a single-particle operator that contains local terms >^(f ) as well as hopping terms 
'fij(t). We note that, because of Eqs. dUi and Q, it follows that 



6 



Michele Fabrizio 



I ^(f)t^(f) %{t) I %(0) = (fo(f) I ^,-(0^^/(0 ^(0 I fo(0) 

= (%(f)|>i(f)l%(f))- (17) 

Seemingly, 

= (fo(f) I ^,-(0^^/(0 ^;(0^^,(0 ^,(0 I %(0> = m) I I fo(0) 

connected • 

The connected term on the right hand side means that we have to extract out 
of ^,(f)^^,(f) a number of fermionic operators, which are to be muhiple of 
two, one of which has to be contracted with yij{t), and the remaining ones with 
^;(f)^^j(f) . By construction, the terms where we extract only two operators and 
average over | 'fb(f)) what remains, will vanish because of Eq. while all the 
others, with four or more operators that are extracted, vanish in the limit of infinite 
coordination number In conclusion, only the disconnect term survives, hence 

(fo(r) I ^(f)t^(f) r,jit) I fo(r)) = mt) I r,jit) I fo(0), (18) 

which, together with Eqs. (T% . imply that 

i{%{t) I ^(f)'^(f) I %{t)) = (fo(f) I ^(f)^^(f) y{t) I %{t)) 

= {%{t) I r(f) I %{t)) = i{%{t) I 'f'o(o>. (19) 

Finally, Eqs. (|4]i and ^ also lead to 

i{%it) I ^(f)^^(f) I fo(f)) = i{%it) I ^,-(f)^^KO I fo(0>- (20) 

f 

As a result, Eq. ( fT6] t can be written as 

^(t) = i{%{t) I %(0) + /I(fo(0 I ^KO^^KO I fo(f)) (21) 



3.1 A more convenient representation 

In order to make it easier the search for the saddle point, it is convenient to fol- 
low the method outlined in Ref. ||2TI . closely connected to the rotationally invariant 
slave-boson formahsm of Ref. 1221 . We assume there exists a set of creation and 
annihilation operators, the natural basis operators c/.^, and d-^, respectively, related 
to the original operators, c^^ and c-^, by a unitary transformation and such that 

{%it)\dld_,\%)=5aiinl{t). (22) 
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We introduce the Fock states in the natural basis 

h-;W> = nfe)""|o), (23) 



such that the matrix P^{t) with elements 



^n}{,n} n ('4(0)"" (1 ^ 5{„}{,„}P°|„j(f), (24) 

a 



is diagonal. We write a generic Gutzwiller projector as 

^iit)= £ ^i^g^ |/;r)(/;{„}|, (25) 

with variational parameters ^',;r{n}(0 that define a matrix <Pi{t), and where | ;;F) 
are basis states in the original representation in terms of the operators ct^. In fact, 
a nice feature of such a mixed original and natural basis representation of the 
Gutzwiller projectors is that one can carry out all calculations without specifying 
what the actual natural basis is: 11211 122ll23l it is just sufficient that this basis exists. 

In this representation, the constraints Eqs. (|4) and (jS) can be simply rewritten 
asED 



Tr(^,(O^^KO^) = 1, (26) 
Tr(«.,(0^^,(0'44) = I dld,^ I %{t)) = Va, (27) 

Tr(4.,(0^^,(0^44) = I did,, I %{t)) =0, Va ^ j3, (28) 



where, from now on, given any operator we shall denote as (9, its representation 
in a basis of states . It turns out that only the constraint (IZTT i requires some care to be 
implemented, while the other two can be implemented once for all at the beginning 
of the calculation]! 



In fact, we can parametrize 



where Ui{t) is a unitary matrix with elements ?7,r{„}, while Pi{t) a positive definite matrix with 
elements -f;{„}{„,} (f). which can be represented as the density matrix of a local normalized state 

lv/(0) = i:^-,w(Oh;W), 

{'-} 

with {Yiit) I 1^,(0) = 1> which automatically fulfills Eq. i26\ . In order to impose the constraint 
{27) it is then sufficient that, for a 5^ j8 

{Mr,{t)\dld.^\Wi(t))=0. 
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In this representation, Eq. Q becomes 

(f (f) I I f (0) = Tx(^Hi)'' OiHt)) , (29) 

hence the average of any local operator can be expressed solely in terms of the 
matrices 4>, without any reference to the Slater determinant. In terms of 4>; one can 
show that 

{%{t) I ^t{t)^Ait) I fo(0) = Tr(^0,(O^^^) • (30) 

Also the effective Hamiltonian can be expressed simply in terms of the 

matrices We define a matrix ^i(f) whose elements are 12211211 



Ru,a{t) = , = Tr 4.,(0^c,„4.,(04 h (31) 

v/«?a(0(l-«l!.(0) ^ ^ 

which, by Eq. dZTl l. can be regarded as functional of 4>,- alone. In terms of those 
parameters. 



N N 



i.j cr,b=l a.p = l 



and we must make sure that this non-interacting Hamiltonian does produces a local 
density matrix diagonal in the dia operators. In conclusion, having introduced the 
matrices <I>i, we can rewrite the Lagrangian (l2ll as 



dt 

+/(%(/) I ^(0) - ('i'oCO I ^* [^(0] I 'i'o(O)- (33) 

We still need to impose the constraint Eq. ( |27] | in a convenient manner In fact, 
what we are going to show now is that we do not need to impose any constraint 
at time f > if that constraint is fulfilled at time f = 0. Since the matrix 4>, is 
variational, we can always write 

with and on the right hand side being independent variables. We assume that 
V,- is a unitary matrix that corresponds to a unitary operator % such that 

'^ya^-^ = L^,a/3^,/J, (34) 
b 



This can be done by regarding | tlie eigenstate of a local many-body Hamiltonian that 

does not contain any term of the form tt^ | /; {«} | t.„ for any | {n}) including the vacumn. 
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It is straightforward to show that 
so that 

where y = H; Therefore the Lagrangian transforms into 
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(35) 

(36) 
(37) 



^(0=E/Tr(4.,'(0t^ 



/Tr(^,'(0t4.,'(0^V,(0 



-Tr(^,'(0^^,-^;(0) 

+i(fo(f) I %{t)) - (%(f) I r jr, [^'(f)] I fo(f)>. (38) 

Since also the Slater determinant is a variational parameter, we can redefine 

\%{t))->y\%{t)), 

where | %'(?)) is still a Slater determinant, because of our definition of Y, and is 
independent of it. It follows that 



dt 



+i{%{t) I %{t)) + i{%{t) I y{t)H{t) I fo'(O) 

-{%{t)\M',[^'{t)\ |fo'(0>, 



(39) 



where the only piece of the Lagrangian that depends explicitly on Y is, being 
unitary. 



5^ 



(40) 



Now, let us assume that 

Yiit) = exp 
It follows that ( l40l ) becomes 



(41) 
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5^ =L<^,a(0 



-Tr 0,'(f)t4.;(r)4,«f,„ + I 4«',a I %'(0> 



(42) 



Since this is the only term that depends on ^,v,, the Euler-Lagrange equation 



d(j)u, dt d^ic 



0, 



impUes that 



dt 
_ d 
~ dt 

- ^ 

~ Tt 



Tr 0;(f)t^;(o4,4 + {%{t) \ dU, \ fo'(O) 



Tr( 0KO'^KO44 ) + ('i'o(f) 1 4^,a I fo(0) 



(43) 



In other words, provided Eq. (I27t is satisfied at f = 0, and Eqs. ( |26] | and (1281 1 are 
enforced by construction, then the constraint (l27t is automatically satisfied by the 
saddle point solution at any time f > . 

In conclusion, under the above assumptions, the only requirement is finding the 
saddle point of the action whose Lagrangian is given in Eq. ( |33] |. Specifically, the 
Slater determinant must satisfy the equation 



i\%{t))=Mi 0{t) \%{t)), 



(44) 



which is just a Schrcedinger equation with a time-dependent Hamiltonian that de- 
pends parametrically on the matrices <Pi{t). These latter in turns must satisfy 



i^^=U,0,{t) + {%{t)' 



Sit) 



dt 



d'p.ity 



%it))=Ht %{t),<P{t) 0i{t), (45) 



which is a non-linear Schrcedinger equation whose Hamiltonian i/, depends not 
only on the Slater determinant | ^o{t)) but also on the same 4',(f) at site / and on 
the 0j(f)'s at the neighboring sites. We note that the time-evolution as set by the 
Eqs. (l44l i and (|45] | is unitary, hence conserves the energy if the Hamiltonian is not 
explicitly time dependent. In other words, one can readily show that 

^ ^ ^(-^(0 I ^ I 'F(O) = j^{%{t) I Mi [0(f)] I fo(0) = 0, (46) 

if I %{t)) satisfies Eq. (l44l l. while ^>, (f ) and ^,(f )^ satisfy Eq. ( |45] | and its hermitean 
conjugate, respectively. If J^{t) is explicitly time-dependent then, under the same 
conditions as before. 
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dE{t) 
dt 



dt 



(f(OI^(f)|f(f))-(%(OI 



tMt) 



dt 



fo(0), (47) 



where the time derivative in the rh.s. only refers to the expHcit time dependence. 
The stationary Hmit of (|44] | and (l45T l, i.e. 



E <P 



A 



dE 4> 



Ur 



30 



0i=H, 



%,0 



(48) 



(49) 



for the lowest eigenvalues E and A corresponds to solving the conventional equi- 
librium problem discussed in section |2] as showed in Ref. Il24l . In particular, the 
Eq. ( |49] | is a self-consistent eigenvalue equation similar to Hartree-Fock, in which 
the Hamiltonian depends parametrically on the same eigenstate that is looked for 

In conclusion, the Eqs. ( |48] ) and (09]) for the stationary condition at equilibrium, 
and the Eqs. (|44] | and (l45T l for the out-of-equilibrium evolution, provide a very sim- 
ple tool for studying the correlations effect in a strongly interacting electron model. 
The method is very flexible; it can deal with many orbitals and also with inhomoge- 
neous situations where the Hamiltonian and/or the initial state are not translationally 
invariant, hence the matrices 0i{t) become site dependent. We stress once more that 
the approach is variational only in the limit of infinite lattice-coordination, otherwise 
it is just a mere approximation without any control parameter, exactly like DMFT 
when it is used in finite and not just in infinite dimensions. 

One aspect worth to be mentioned is that within the Gutzwiller approximation 
two different types of dynamical degrees of freedom seem to emerge. One is pro- 
vided by the Slater determinant with its evolution (|44] |. It is commonly believed that 
this set just describes the quasiparticle degrees of freedom. In addition, the matri- 
ces <Pi introduce other local degrees of freedom with their own dynamics set by 
Eq. i45[ . It is tempting to associate them with the incoherent excitations that coexist 
with the coherent quasiparticles in the presence of interaction, and which become 
the Hubbard bands near a Mott transition. lfT6l Within the Gutzwiller approximation, 
coherent and incoherent excitations are coupled to each other in a mean field like 
fashion, which provides a very intuitive picture although it misses important dissi- 
pative mechanisms. In what follows, we shall provide additional evidences that <Pi 
are indeed related to the Hubbard bands. 



3.2 A simple case study 

Before concluding this section, we think it is worth showing how the equation sim- 
phfy in the frequent and relevant cases in which the point symmetry of the Hamil- 
tonian already determines the local orbitals in which representation the local single- 
particle density matrix is diagonal, i.e. 
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(f(f)l4,c.J'f'(0) = 5„,«,„(0. (50) 
In this case, where natural and original basis coincide, hence also 

{m)\cic.,\%it)) = d„„nlit), (51) 
the expression (ISTT l further simplifies into 



Rl,it) = dab I Tr(^,(0^4^'(0q„) ^ RUt) Sat- (52) 

^«?«(0(i-«!i(0) 

Because of the constraint Eq. ( |27] i. we can equivalently regard 

ni{t)^Ti(4>i{t)'<P,{t)nu), (53) 



as functional of ^, (f ), rather than of the Slater determinant, hence it follows that 

dRUt) _ 1 



2«?„-l 



+RUt) w '" X Mt)n,a, (54) 



2«?„-l 



+^'"(0 . / n \ (55) 
2«L(l-4-(0) 

If we consider the Hamiltonian 
then 



= I Y.'tj{RiMRjb{t)cic,„+H-c) , (57) 

'7 ab 

so that, through ( |33] |. the Slater determinant satisfies that Schroedinger equatioij^ 

/|'f'o(0)-^^*(0|foW)- (58) 



^ Once again, we must make sure that the effective Hamiltonian ^(f), Eq. l |57K is such that the 
local density matrix remains indeed diagonal in the operators ^ . 
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If we define 

A^t) = E t^^Rjbit) (fo(f) I c^,,cj, I %{t)), (59) 

jb 

then 0i satisfies the matricial Schroedinger equation 



I- 



" /^°(0(i-<(0 



2ni - 1 



+ E /?*,(OA«(f) + c.c. 4.K0«,«. (60) 

« ^ ^2««(l-n«(0) 

The equation for can be obtained simply by the hermitean conjugate of ( |60] l. We 
can readily demonstrate, through (|52t . that 



" y«?„(o(i-<(o) 
y«?«(o(i-«°(o 



: Tr 0, 



^ : Tr( m^clmc^t 

^«i;(o(i-<(o) 



y„0,(r)(l-„0(0) 

= /J,-i(0*Ay,(r)-/?,-,(r)Afc(0* 

= Ef^"(/?,7.(r)*/?,„(f)('f{,(0 I 4c^„ I fo(0> -cc.) 
if 

= (%(f) I [«,fo,J^*(f)] |fo(f)) 

= i^{%{t)\n,i,\%it)), (61) 

which explicitly proves that the constraint is indeed conserved by the above dynam- 
ical evolution. 
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4 Quantum quenches in the half-filled Hubbard model 



Armed with all previous results, we can start investigating the simplest possible out- 
of-equilibrium evolution in the single-band Hubbard model at half-filling. For sake 
of simplicity we shall ignore magnetism, hence assume spin SU (2) invariant | ^Q{t)) 
and 4>, . In this case, natural and original bases coincide, hence we can use the results 
of section [J!2l We choose as a local basis that of an empty site, | 0), doubly-occupied 
site, I 2), and singly occupied site with spin up, |t), or down, \\.). We take for <Pi with 
elements <Pirr' with r,F' = 0,2,t,| the SU (2) and particle-hole invariant form 



1>, = ^ 
V2 



(^m \ 

4>,22 


V 0,;;/ 



(62) 



with 0,00 = 't'iii = ^io and = = ^>,i . All constraints Eqs. (|26ll-(l28]l, with 
n^-^{t) = n'^^{t) = 1/2 Vf, are satisfied provided 



|^,-or + |^n 1=1. 
With the above parametrization the Eq. (|52] | becomes 

Ri^it)* ^ Ri<^{ty = Riity - 0,-o(r)*0,i(O + ^n(O*^,o(O e Me. 
Given the original Hamiltonian 



then 



Mit) = ^ t,j [Ri{t)Rj{t)clcj,+H.c.), 



(63) 



(64) 



(65) 



(66) 



and the Slater determinant is the solution of the Schroedinger equation ( |58] l. In this 
case in which a =t,4 and spin symmetry is preserved, the parameter defined in 
Eq. ^ 

At(0 = Hit) = ^ = f„^;(0 (fo(0 I I fo(0> e Re, (67) 

is real. Therefore the equation of motion (|60) becomes 



i^n{t)=2Ai{t)<P,a{t). 



(68) 
(69) 



We note that if we imagine the spin- 1/2 wave-function 
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\^,{t))^^n{t)m+^io{t)m (70) 
solution of the Schroedinger equation of the spin Hamihonian 

^*isme = Y.^[l-cyf)+2A,it)a[, (71) 

that describes independent spins in a uniform magnetic field —U/4 along z and a 
site and time dependent field 2zi,(f) along x, we would get exactly the equations 
dMll and with 

Riit) = {<P,\a'\<Pi), (72) 

implying that the field 2Aj{t) is self-consistently determined by the same spins. This 
observation is not a coincidence, as we shall discuss later 

Before analyzing a simple case of out-of-equilibrium evolution, let us consider 
the stationary limit, which, as we discussed, defines the equilibrium conditions. In 
this case it is likely that the lowest energy state is homogeneous, namely invariant 
under translations, hence Ri = R, V/. The stationary solution of Eq. (ISST i is just the 
ground state of the hopping Hamiltonian with renormalized hopping parameters 
t^,jj = R^tij and energy per site R^eo < 0. Therefore 4,- = A = Rsa, for all /, hence 
the Eqs. ( |68] l and ( |69] l in the stationary limit become simply (we drop the site index 
as all sites are equivalent) 

A0o = ^<Po + 2eoR0u (73) 
= 2eoR<pQ. (74) 

We write <Po = sin 9 /2 and <Pi = cos 9 /2, so that the wavefunction is normalized, 
hence R = sm9. The eigenvalue problem is solved if 

for f/ < J/c = 8|eo|, in which case 

U U , , 

A = - + 2eo = --2|eo|, 

otherwise, for U > Uc, the solution is = with energy A = 0. Indeed, Uc can 
be identified as the critical repulsion for the Mott transition within the Gutzwiller 
approximation, because, for U > Uc , R = sin 0=0, hence the hopping energy van- 
ishes. We observe that the highest energy eigenvalue at self-consistency is 

A' = ^ + 2|eo|, 

for U < Uc, and A' — U /2 above, resembling much what we would expect for the 
location of the Hubbard bands. 
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Let us come back to the out-of-equilibrium evolution, and suppose we start at r = 
from the ground state of the non-interacting Hamiltonian, which is just the ground 
state average of the hopping with energy per site Co < introduced above, and total 
energy £0 < 0. This corresponds to assuming that 4',o(0) = ^",1 (0) — 1 /a/2, hence 
Ri{0) = 1, V;, and | ^(O)) being the uniform non-interacting Fermi sea. Since trans- 
lational symmetry remains unbroken during the time evolution, Ri{t) V/ and 

Vf > 0. It follows that ^ (f ) remains the same tight-binding Hamiltonian as at f = 0, 
just renormalized by the overall factor R{t)^. As a result, the Slater determinant evo- 
lution is trivial, 

I %{t)) = e-'^o/o I «fo(0)), (76) 



hence zi,(f) ^ R{t)eo. Therefore the equations (|68] | and (|69] | become for any site 
equal to 

'^o(0 - ^Mt) + 2eQRit)<Pi{t), (77) 
i<t>i{t)^2eoRit)0oit), (78) 
with/?(f) = <Pi (f )*^>o(f) +C.C.. If we set 

{<P{t) I a' I <P{t)) = sin0(f)cos^, (79) 

(0(f)|(T-^'|0(r)) =sin0(Osin^, (80) 
(0(r)|(j-^|0(r)) =cos0(r), (81) 

then, through Eqs. ( |77] i and ( |78] l. we find the following equation of motion for (j){t): 



0(f) = ±yf/2_i6e2sin2 0(f), (82) 

which is just the equation of a pendulum. In particular, if U < 4|£o|, 0(f) oscillates 
between ±0Max, where 

0Max = Sm 



4|eo| 

On the contrary, when U > 4|eo|, 0(f) increases indefinitely. In other words, the 
quench dynamics displays a dynamical critical point at [/* = 4|£o|.|9] We observe 
that U^: is just one half of the critical Uc that we found previously at the Mott transi- 
tion within the Gutzwiller approximation. Remarkably, an abrupt change of dynam- 
ical behavior near Uc/2 has been observed also in Ref. [7| within a time-dependent 
DMFT simulation of the same quantum quench as above. Given the very crude ap- 
proximation in using a Gutzwiller wavefunction with respect to the exactness of 
DMFT in infinite coordination lattices, such an agreement is indeed quite remark- 
able. 
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In sectionOwe already noticed that the variational degrees of freedom introduced by 
the projectors are promoted to the rank of true dynamical degrees of freedom in 
the time dependent extension of the Gutzwiller approximation. Moreover, in section 
|4]we found that in the simple case of a single-band Hubbard model at half-filling, 
these new dynamical objects resemble spins in a self-consistent magnetic field, see 
Eq. ( fTTT i. In what follows we will put such an analogy on a more solid basis, although 
the demonstration applies rigorously only to few simple cases. The outcome will 
be a theory that looks similar to the so-called slave-spin representation recently 
introduced lfT2l [T3] [T4l [TSl as an alternative approach to slave-boson theory. 



5.1 SU {N) Hubbard model at half-filling 

We note that, at given | the Gutzwiller wave-function | f ) in Eq. (|2]i actually 
defines a whole set of wave-functions, each identified by the projectors that act 
on each site ;. Let us assume there exist a whole set of projectors ^,„, that satisfy 



{% I Kr^.n^L^.ta' I '^'o) = d„,„ {% \ c^c,,^ | %) , ya,b and V(T,(t', (84) 



where we distinguish between orbital indices, a,b ~ 1 ,A^, and spin indices, cj 

and a'. It is straightforward realizing that these conditions allow to evaluate, along 
the same lines previously outlined, also matrix elements between different wave- 
functions. In this way, one can get the matrix representation of the Hamiltonian on 
such a subspace of wave-functions, whose diagonalization provides not only a better 
estimate of the ground state energy but also gives access to excited states. 

The Hamiltonian we shall consider is given by ([TJ with diagonal nearest neighbor 
hopping -8aht/^/z and 



where «, = LaaSaCT'-iaCT' density corresponds to electrons per site, i.e. 

half-filling. The model therefore is invariant not only under spin SU (2) but also 
orbital SU{N), in fact it is invariant under the large U{2N) symmetry group. We 
shall therefore assume that the wave functions | f ) and | %) are invariant under 
such a large symmetry. We define =S„, the projection operator at site / onto states 
with n electrons. If we choose as local basis the Pock states | /;{«}) identified by 
the occupation numbers niaa = 0, 1 in each orbital and spin, i.e. 



(fo I ^,;,^,„ I %) = 5, 



(83) 




(85) 
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then 

= ^ ( " ^ ) I {«}> (''; {«} I ■ 

From the invariance properties of the Slater determinant | ^q) it follows that 



{% I I fo) = i P,f ) ^ Pi'') = (86) 



as well as that 

I ^,-„ I %^ 

4' 

where fi^^ is the distribution probability of the local occupation number on the 
uncorrected wavefunction. The most general Gutzwiller projector satisfying (|4]i and 
(|5]l can be written as 



2N 



= ^ ^ (87) 

where 



=0 Jp(0) 



AT 



s=-N 

and I ^lif 1 = 1 01 -s I - In fact, we can regard as the wavefunction components of 
fictitious spins of magnitude S ^N, one at each at site /, 

.v=-S 

which we shall intentionally denote as slave spins as they are closely related to the 
slave-spin representations of Hubbard-like models. ifTZl [131 [T4l [TSll 

The renormalization factor defined by Eq. (|9) is in this case diagonal, R-iat — 
Ri 5ab, and simply given by 



s-i 



1 



= ^{0^\S+\0i). (88) 

More generally, the matrix element of the fermionic creation operator cj^^ between 
two wave-functions, | f ) and | f'), with local projectors and at site hence 
slave spin wave functions | ^>,) and | 0'-), respectively, has the very transparent 
expression 



s 

Seemingly, the matrix element of the local repulsion reads 
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I I («,■ - A^) V> = I (^H (S'f I , (90) 
where is the z-component of the slave spin operator S. In conclusion, we find that 

{W\J^\y")=--^ ^ ({<P,\S+\<P[){<Pj\S-\<P'^){%\cl^c^^J%)+H 

indeed a very suggestive result. Notice, however, that the slave spin wave-functions 
are not completely free, because they must correspond to Gutzwiller projectors sat- 
isfying ^ and (EDi. 



5.2 Slave-spin basis 

Therefore, to make Eq. (l9ll suitable for calculations, we still need to identify a 
proper set of Gutzwiller projectors satisfying Eqs. (|8J| . (l84l i. A possible choice is 




with m < A^. In principle we could have also chosen the combination ( |93] l with the 
minus sign instead of the plus, but not both as they are not orthogonal in the sense 
of Eq. ( |84| |. In other words, not the whole slave-spin Hilbert space is allowed, but 
only a subspace | (m)), with m = 0, ... ,5: 

l(0)>^|0), (94) 
|(m>0)) = -J=(|m)+|-m)), (95) 

which we shall denote as the physical subspace, still to keep contact with the jargon 
of slave-boson theories. 

We note that the action of the raising operator projected onto the physical 
subspace, i.e. 

S+l(0))-\/S£+II|,i),, (96, 

s-i(i))^/sk|,„»+;s£MH|p,), 
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r,+ I / -,Nx I S(S+ 1) -m(m+ 1) , , 
5+|(m>l))~ K'w+l)) 



|(,.-1)>, (98) 



4 

is just the same as the action of 5* without any restriction. Indeed 



I (o)> ^ 5-^- 1 0) = ( 1 1>+ 1 -i>) = yfcll I (D), 

5^ I (1)> = ^ (v/5(5+ 1) - 2 I 2) + v/5(5TT) I 0) 



+ VS{S+l) I 0) + ^S{S+l)-2 I -2) 



and also ( |98l ) follows trivially. Actually, the physical subspace is invariant under the 
action of 5^, therefore, using the latter instead of 5+, we are allowed to release the 
constraint and work in the full Hilbert space of the slave spins, since we expect the 
ground state to contain | 0), hence to occur within the physical subspace of Eqs. ( |94l i 
anddlSl). 

In conclusion, we can rewrite (l9Tl i as 

<i.j>aa 

+ yL(^H(5-f 1^,'), (99) 

without any condition to be imposed on the slave spin wave-functions. We finally 
note that Eq. ( |99] l is just a matrix element of the Hamiltonian 

= E SlS^UaCjao+H.c)+jj:iSif, (100) 

<ij>aa ^ ^ ^ i 

which describes electrons coupled to slave spins of magnitude 5 = A^. In this repre- 
sentation the slave spins are not subject to any constraint. 

We note that the Hamiltonian dlOOl ) resembles much the slave-rotor representa- 
tion for the multi-orbital Hubbard model of Ref. ||25l , with however a major differ- 
ence. In fact the Hamiltonian dlOOl ) possesses only a discrete Z2 gauge symmetry, 
unlike the slave-rotor Hamiltonian that has a larger U{1) gauge symmetry. This dif- 
ference has some important consequences that we discuss below. 
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5.3 The Mott transition 

The great advantage of the representation (1100) is to make the Mott transition ac- 
cessible akeady within the mean field approximation. The simplest mean-field ap- 
proach amounts to assume a factorized variational wave-function | f ) = | electrons) x | 
slave-spins) . The minimum energy is obtained by choosing | %) the Fermi sea of a 
simple tight-biding Hamiltonian. If we define 

' V ^ <i.j>aa 

the hopping energy per site of the state | %), then the slave-spin wavefunction must 
be the ground state of the Hamiltonian 

^smg = -4 - L 5?^)+ T L (S^)' ■ (101) 

This spin Hamiltonian has a discrete Z2 symmetry S\ — 5/, V/, which is spon- 
taneously broken at small U /J, i.e. {Sfj is non-zero and corresponds to the order 
parameter, and restored only above a quantum critical point. This Ising-like transi- 
tion corresponds to the Mott transition in the original interacting model. In fact, the 
physical electron translates in the model ( II 00b into the composite operator S]c\^ 
hence, within mean-field, the long distance density matrix 

lim (4^,,)^ lim (5f45)c.J = lim (5?5j) (4c^.,). 

The average over the electron wave function, which is the ground state of the hop- 
ping, is long ranged. Therefore the long distance behavior of the physical electron 
density matrix depends critically on the slave-spin correlation function. In the sym- 
metry broken phase, 

lim {S',S') {S-y ^ 0, 

\i-j\^oo 

hence the physical electron density matrix is long ranged, as we expect in a metallic 
phase. On the contrary, when the symmetry is restored, then {SJS'^) vanishes ex- 
ponentially for \i — j\ °°, transferring such an exponential decay to the physical 
electron density matrix, which therefore does not describe anymore a metal phase 
but rather a Mott insulating one. It is important to notice that, in the actual slave- 
spin model (II 00b . a finite order parameter {S'J) corresponds to a phase with broken 
Z2 gauge symmetry, which is possible in spite of the Elitzur's theorem li26l because 
we are working in the limit of infinite lattice coordination. ||27l[28l We also observe 
that in the symmetry broken phase there are not Goldstone modes because the sym- 
metry is discrete, unlike what predicted by the slave-rotor mean field theory, Ii25l 
where these gapless modes are expected and associated with the zero-sound. 

The location of the Ising critical point of the slave-spin Hamiltonian dlOOb can be 
determined approximately by assuming that it occurs for large enough [/'s so that it 
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is safe to keep only states with 5-" = 0,±1. We denote 

lt) = |0), 




as the two states of an Ising variable, and introduce Pauli matrices in this subspace. 
We find that the operator 5'^' in this subspace acts like y/S{S+ l)/2, while {S~)^ 
like (1 — C7^') /2, so that (llOlb can be rewritten as 

^3,ng - -/ \ + ^ E (1 - ^1) , (102) 

namely like a simple Ising model in a transverse field. We note that for the case 
5=1, the Hamiltonian jlQH coincides with the slave-spin representation of the 
single-band Hubbard model. |[T2l |l4l |T5l [TOl 

The model (|102| i has indeed a quantum phase transition that separates a ferro- 
magnetic phase, {a'D ^ 0, for small U jJ, from a paramagnetic one, (aj) — 0, for 
large JJ /J. This transition is actually the Mott transition in the slave-spin language 
and, within mean-field, it would occur at a critical 

(103) 

We note that in the single-band case, S ~ \,Uc coincides with the value obtained 
previously. 

Apart from making the Mott transition accessible by mean-field, the effective 
slave-spin model also uncover new dynamical excitations that it is natural to asso- 
ciate with the Hubbard bands. Indeed, the models dlOlb and its simplified version 
(IIO2I 1 display a spin-wave branch that becomes soft only at the transition. For very 
large U , the excitation energy becomes ^ U /2, just the location of the Hubbard 
bands. Needless to say, the mean field dynamics of these spins corresponds to the 
dynamics of the matrices 4>, that we introduced previously. 



5.4 Away from half-filling 

We can repeat all the above calculations even away from half-filling. In this case, the 
slave-spin wave-functions | ^a) in the physical subspace must satisfy the conditions 

(0a I ^p) = 5„/3, (^a I S'' I f^) = (104) 

where 5 =n — N is the doping away from half-filling. The expression of the varia- 
tional energy is modified into 
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(5- - 5 ) Vi <,-,;>a« V ' ' / 

+ yL(«'-l(5f 1^/). (105) 

As before we need to identify the physical subspace for the slave-spins. 

The simplest case is when the average occupancy n is integer, hence 5 is integer, 
too, which requires more than a single band, i.e. N > I. Let us further assume U 
large, so that we can just focus on the two physical states 

\i)^^{\d+i)+\d-i)), 

which corresponds to the assumption that a kind of particle-hole symmetry is recov- 
ered close to the Mott transition. The raising operator projected onto this subspace 
has the action 







-1) 


2 




-l)-5{5- 


-1) 



^"i;)^V 2 it)-(«-^) It), 

hence 5+ ~ ac* — /)3(7^, with a >| j3 |. It follows that the Ising variables are de- 
scribed by the effective Hamiltonian 



with / being the average hopping per site of the Fermi sea with average occupation 
n. This model still has a phase transition between a ferromagnetic phase with (c'^) ^ 
and a paramagnetic one. Within mean field, the critical interaction strength is now 

Uo^^Js^' ^^^^^ 

and is shifted to lower values of the interaction as 8 increases. Once again, the spin- 
wave spectrum of the Ising model (IIO6I 1 can be interpreted as the spectrum of the 
Hubbard bands. 

If the filling is not an integer or the enlarged SU{2N) symmetry is lowered, the 
above construction does not work anymore because we cannot define in general 
more than a single Gutzwiller projector satisfying both dSST l and ( l84l i. In other words, 
while for integer fillings and SU{2N) symmetry we can associate the dynamical 
variables <Pf with auxiliary spin operators, which allows for instance to improve the 
Gutzwiller approximation by including systematically quantum fluctuations, away 



24 



Michele Fabrizio 



from such a high-symmetry points we are unable to make such a simple identifica- 
tion, hence we must limit our analysis to the mean field dynamics of <Pi. 



6 Conclusions 

In this paper we have shown in detail how one can access by simple means the out- 
of-equilibrium time evolution of a Gutzwiller-type variational wave function. The 
approach is rigorously variational in the limit of large coordination numbers, other- 
wise can be regarded as the dynamical counterpart of the widely adopted Gutzwiller 
approximation. The method is really simple to implement and very flexible. It is 
apt to cope with weak non-equilibrium compatible with linear response, but also 
with strong out-of-equilibrium conditions like sudden quantum quenches. It can de- 
scribe single- and multi-band systems, as well as homogeneous and inhomogeneous 
models. 

The key feature that distinguishes the present method from the conventional time- 
dependent Hartree-Fock is the emergence of two distinct types of excitations that 
control the time-evolution of the wave function. One corresponds to the particle- 
hole excitations of the guiding Slater determinant, just like in the time-dependent 
Hartree-Fock, and is supposed to describe coherent quasiparticles. In addition, new 
local dynamical degrees of freedom emerge, which can be associated with the Hub- 
bard bands and that are promoted to the rank of genuine excitations with their own 
dynamics. Within the Gutzwiller approach the Hubbard bands and the quasiparti- 
cles are mutually coupled in a mean-field like fashion, i.e. each of them generates 
a time-dependent field that acts on the other In spite of such an approximation, 
the dynamical behavior that follows is quite richer than in Hartree-Fock. We have 
shown just an example of such a richness, namely the dynamical transition that oc- 
curs in the single-band Hubbard model at half-filling after a sudden increase of the 
repulsion. |j9l 

Finally, we have shown that it is possible to extend the variational approach 
to a multi-configurational wave function that comprises a linear combination of 
orthogonal Gutzwiller-type of wave functions. Such a multi-configurational vari- 
ational method can be worked out analytically only in specific cases, specifically 
for integer fillings. Nevertheless it is quite instructive since it demonstrates that 
the above discussed time-dependent Gutzwiller approach is nothing but the mean- 
field approximation applied to the actual Hamiltonian dynamics within that sub- 
space of orthogonal Gutzwiller wave functions. Remarkably, the Hamiltonian pro- 
jected in that subspace resembles the slave-spin representations of correlated elec- 
tron models. lfT2l [131 [141 [TSl thus providing a very intuitive picture of these theories. 
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